# include <stdio.h>
# include <stdlib.h>
# include <math.h>

/* cmtrx[nxpt-1][nstn].tmtrx[nxpt-1].x[nxpt].xpos[nxpt-1]. */
static short cmtrx_interp(unsigned short nxpt, double *x, double *xpos,\
 unsigned short nstn, double **cmtrx, unsigned short *tmtrx);

int main()
 {
  short nxpt = 14;
  short nstn = 5, i,j;
  double x[15], u[15], val, xpos[14];
  double **cmtrx;
  unsigned short tmtrx[14];

  cmtrx = (double **)malloc(sizeof(double *)*(nxpt-1));
  for(i = 0; i < nxpt-1; i++) cmtrx[i] = (double *)malloc(sizeof(double)*nstn);

  for(i = 0; i < nxpt; i++)
   {
    x[i] = asin(-0.91*cos(2*acos(0)*i/(double)nxpt))/asin(0.91);
    u[i] = exp(-x[i]);
    if(i > 0) xpos[i-1] = /*x[i-1] - */0.5*(x[i-1]+x[i]);
   }

  cmtrx_interp(nxpt,x,xpos,nstn,cmtrx,tmtrx);

  for(i = 0; i < nxpt-1; i++)
    {
     for(val = 0, j = 0; j < nstn; j++) val += cmtrx[i][j]*u[tmtrx[i]+j];
     printf("%0.15E : u[n=%d]\n%0.15E : Real Value\n\n", val, i, exp(-xpos[i]));
    }

//   for(i = 0; i < nxpt-1; i++) { for(j = 0; j < nstn; j++) printf("%E\t", cmtrx[i][j]); printf("\n"); }

  for(i = 0; i < nxpt-1; i++) free(cmtrx[i]); free(cmtrx);

  return 1;
 }


static short cmtrx_interp(unsigned short nxpt, double *x, double *xpos,\
 unsigned short nstn, double **cmtrx, unsigned short *tmtrx)
  {
/* cmtrx[nxpt-1][nstn].tmtrx[nxpt-1].x[nxpt].xpos[nxpt-1]. */
   unsigned short irow, i, j;

   if(xpos == (double *)0) return -1;
   if(nstn >= nxpt) return -2;

   for(i = 0; i <= nxpt-2; i++)
       {
     /* Marking begining of the stencil on the global co-ordinate. */
        if(i <= (nstn/2)-1) tmtrx[i] = 0; /* Zone-1 */
        else if(i <= nxpt-(nstn/2)-2) tmtrx[i] = i-(nstn/2)+1; /* Zone-2 */
        else tmtrx[i] = nxpt-nstn; /* Zone-3 */
       }

   for(irow = 0; irow <= nxpt-2; irow++)
   for(j = 0; j < nstn; j++)
   for(cmtrx[irow][j] = 1.0, i = 0; i < nstn; i++) if(i != j)
      cmtrx[irow][j] *= (xpos[irow]-x[tmtrx[irow]+i])/(x[tmtrx[irow]+j]-x[tmtrx[irow]+i]);

   return 1;
  }
